Spin-supersolid phase in Heisenberg chains: a characterization via Matrix Product 

States with periodic boundary conditions 
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By means of a variational calculation using Matrix Product States with periodic boundary con- 
ditions, we accurately determine the extension of the spin-supersolid phase predicted to exist in 
the spin-1 anisotropic Heisenberg chain. We compute both the structure factor and the superfluid 
stiffness, and extract the critical exponents of the supersolid-to-solid phase transition. 
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A phase of matter where diagonal (solid) and off- 
diagonal (superfluid) long-range order coexist is named 
supersolid. Since its original prediction,— the search for 
this phase has attracted the attention of a growing num- 
ber of experimental and theoretical physicists^ However, 
despite this great effort, the supersolid phase has, to date, 
eluded a firm experimental confirmation. This is due to 
the fact that the stabilization of such a phase arises from 
the combined action of two mutually exclusive effects: on 
one side, the solid order requires a well defined spatial ar- 
rangement of the atoms in real space; on the other side, 
the superfluid order requires the atoms to be delocalized 
and condensed in a macroscopic quantum state. 

The first, and probably most prominent, candidate for 
the experimental realization of a supersolid phase is He»2 
More recently, the trapping of cold atoms in optical lat- 
tices has stimulated the search for such exotic phase in 
these systems (see Refs. [J and references therein). 

Furthermore, in strict analogy with what postulated 
in the fields of quantum fluids and cold atomic gases, a 
spin-supersolid phase can be defined also in the context 
of quantum magnets, in association with a simultaneous 
ordering along the z-direction at finite momentum and of 
a breaking of U(l) symmetry in the xy-plane. Examples 
of such phases have been found^— in S — 1/2 spin-dimer 
model on a square lattice, where extra singlets delocalize 
in a solid background via correlated hoppings^ and in 
5 = 1 systems 

The spin-1 Heisenberg chain with a single-site uniax- 
ial anisotropy in a transverse magnetic field is what we 
study in the present paper. For this model, Sengupta 
and Batista^ predicted a spin-supersolid phase for in- 
termediate values of the external field and of the uni- 
axial anisotropy. Their analysis of the phase diagram 
was based on the derivation of an effective model and 
on Quantum Monte Carlo simulations. Further con- 
firmation using Density Matrix Renormalization Group 
(DMRG) was reported in Ref. In these last works, 
the existence of the supersolid phase was inferred by an 
analysis of the magnetization profiles. However, due to 
the intrinsic limitation of standard DMRG techniques to 
the case of Open Boundary Conditions (OBC), it was 
impossible to access the superfluid order parameter with 
such kind of algorithm. A detailed numerical analysis of 
the supersolid phase would indeed require the simulta- 



neous study of both diagonal and off-diagonal orderings. 
The Matrix Product States (MPS) approach to DMRG^i 
with its recent generalization to study efficiently one- 
dimensional systems with Periodic Boundary Conditions 
(PBC)fi^— appears to be an ideal tool to determine such 
parameters. Here we exploit this fact to address both 
the diagonal and off-diagonal order parameters for the 
spin-1 Heisenberg model of Refs. l9lflu this allows us to 
directly access the so called spin stiffness of the system, 
and therefore to accurately locate the supersolid phase. 

The model under investigation is governed by the fol- 
lowing spin-1 Heisenberg Hamiltonian 
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where (with a = x,y,z) are the spin-1 operators for 
the j-th site, while are the associated raising/lowering 
operators; PBC are imposed by requiring S^ +1 = S" 
(N is the number of sites in the chain). Notice that, 
in addition to the exchange coupling J and the mag- 
netic anisotropy A, the model also includes a coupling to 
an external transverse field B and a single-site uniaxial 
anisotropy of strength D. Hereafter we set J = 1, thus 
fixing the energy scale. Furthermore, following Ref. 
we set A = 2D. Units of % = kb = 1 are used. 

The phase diagram described by the model in Eq. ([1]) 
is quite rich (see, e.g., Fig. [T]). For large values of the 
anisotropy A, it goes into a spin-gapped Ising-like phase 
showing long-range diagonal order. On increasing the 
external field, there is a transition to a superfluid phase 
characterized by a finite spin-stiffness. At larger values of 
B, the system goes into a fully polarized state (not shown 
in Fig.[l|. In between the spin-gapped and the superfluid 
phase, a spin-supersolid was shown to exist^ possessing 
simultaneously diagonal and off-diagonal ordering. We 
concentrate on this specific configuration. 

The solid ordering can be detected by an analysis of 
the spin-structure factor, defined as 
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Figure 1: (color online). Extension of the supersolid phase 
(cyan region) in the A — B plane, for a one-dimensional 
anisotropic spin-1 Heisenberg Hamiltonian in a transverse 
magnetic field, and single-site uniaxial anisotropy defined in 
Eq. |TJ). The value of uniaxial anisotropy is fixed at D = A/2. 
The phase boundaries are located by evaluating the region of 
parameters for which the solid order parameter Osdw and the 
spin-stiffness p s of Eq. <(3j were simultaneously different from 
zero. For A > 6, where the transition is between the solid 
and the supersolid, the vanishing of the superfluid stiffness 
is an excellent indicator of the supersolid boundaries. For 
smaller values of A the transition is to a superfluid phase, 
therefore here the boundary of the supersolid is determined 
by the vanishing of the solid order (blue squares), while the 
spin stiffness vanishes at larger values of the external field B 
at the boundary between the superfluid and the spin-gapped 
phase (open circles). The two dashed lines are the result of 
effective low-energy models, valid for A ^> 1.— The dotted 
lines are directly taken from the simulations of Ref. [9| and 
separate the solid phase from the superfluid region at large 
values of B and small values of A. In the figure S = Solid, 
SS — Supersolid, SF = Superfluid. 



at momentum q 



A solid order parameter can be 



defined as Osdw — limjv- 



N 



indeed non zero 



values of this quantity indicate that the dominant cor- 
relations have a Spin Density Wave (SDW) character. 
Off-diagonal order instead is detected by the superfluid 
stiffness, denned as 
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where Eq(4>) is t ne ground state energy of the chain with 
twisted boundary conditions, or equivalent!}*^ in the case 
where J — > Je % ^l N ' . For PBC p s quantifies the system's 
response to an infinitesimal magnetic flux which is added 
through the ring. Vice-versa for OBC it nullifies, since 
the twist <f> can be wiped out by a gauge transformation. 
The simultaneous nonzero values of Eqs. © and ([3]) sig- 
nal the supersolid phase. Our investigation leads to the 




Figure 2: (color online). Spin stiffness for model |T} with 
A = 6 and D = A/2 = 3, in a parameter range where the 
system is in a supersolid phase. Parameters used in the MPS 
variational wavefunction for the various sets of data are as 
follows: for m = 10 we performed n s = 30 sweeps, with 
truncation parameter p = 25, 20, 15, respectively for TV 
60. 90, 150; for m = 15, TV = 150 we used p = 30, n s = 40; for 
m — 20, 40 we respectively took p — 45, 60 with n s = 50. We 
kept s — 40 in all the cases except for m = 40, where s = 60, 
obtaining comparable precisions in the energy fluctuations for 
each of those parameter settings. 



result summarized in Fig. [TJ In the following we provide 
detailed evidence of this result. 

Our algorithm is based on Refs Hlilil . where details of 
the implementation can be found. We considered chains 
up to TV = 180, where no finite-size effects could be de- 
tectable for our precisions. The dimension of the matri- 
ces used in the MPS ansatz with PBC was taken up to 
to = 40, while the minimization of the ground state en- 
ergy was obtained by optimizing the structure site by site, 
sweeping through the ring in a circular fashion with a suf- 
ficient amount n s of sweeps. As discussed by Pippan et 
al.r^ an important speedup in the code can be achieved 
by introducing a factorization procedure for long prod- 
ucts of MPS matrices, which reduces the computational 
effort. Intuitively this is justified by the fact that, for 
large chains, the local physics of the system is not af- 
fected by the properties of the boundaries. The degree of 
this factorization is characterized by two truncation pa- 
rameters p and Sj 1 ® that in our simulations were taken to 
be 10 < p, s < 60 (for a formal definition of these quan- 
tities we refer the reader to Ref. Il6|). We checked that 
our choice of m and p would guarantee the convergence 
of our results. 

For the calculation of the stiffness, we computed the 
dependence of the ground state energy as a function of 
the twist and then fitted the curve with a quadratic law 
Eo(4>) = Eq(0) + C20 2 , obtaining the prefactor c-i which 
is directly related to the stiffness: p s = 2Nc2- The deter- 
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Figure 3: (color online). Scaling of the spin stiffness for the 
supersolid-to-solid transition at A = 6. The critical point is 
estimated to be B c = 8.5052 ± 0.0005. Black circles is the 
same data set for m — 10, N = 90 in Fig. [2] Red squares 
are for m = 15, N = 180, with p = 30, s = 40, and n a — 
50. The scaling is compatible with a power-law behavior of 
exponent j3 s = 0.5 (dashed blue line, plotted as a guidance). 
The power-law fits of the two data sets until the vertical blue 
line, respectively giving /3 S — 0.521 and 0.511, confirm this 
prediction. 
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Figure 4: (color online). Spin structure factor S zz (ir) as a 
function of the system size N for A = 6, D = A/2, while 
B — 8.8 (left panel) and B = 15 (right panel); this has been 
obtained with MPS variational technique with OBC. Data 
are rescaled over N. Dashed lines are linear fits of the three 
points at the largest sizes, for m = 100. A finite value of 
Osdw ~ 0.046 ± 0.001 can be obtained by extrapolating the 
N — > oo value in the left panel. On the other hand, in the right 
panel a value of Osdw ~ 1-18 x 10 -4 ± 10 -4 is extrapolated 
at the thermodynamic limit. This corresponds to a vanishing 
solid order parameter, within numerical accuracy given by the 
linear fits. 



mination of the boundary for the solid order turned out 
to be more demanding, due to the necessity to measure 
long-range spin correlations, i.e., the quantities (SjSj +r ) 
of Eq. @, for r ^> 1. Contrary to the evaluation of 
ground-state energies that enter in Eq. (J3J , this generally 
requires a high degree of accuracy in the MPS represen- 
tation of the ground state, thus implying large values of 
to. To enhance the precision, we hence used the fact that 
the solid order in the bulk of the system is not qualita- 
tively affected by the choice of OBC or PBC, and ran 
simulations using MPS with OBC 11 which allows one to 
work with matrices of larger dimension (i.e., with m of 
order 100). We also carefully checked that the obtained 
results were not plagued by finite-size corrections. 

Our findings are summarized in Fig. [TJ which details 
the phase diagram of the system obtained by computing 
the solid and superfluid parameters Osdw and p s for 
different values of B and A. Consider first the results 
we obtained for the superfluid stiffness focusing on a sin- 
gle value of the anisotropy, say A = 6. The behavior of 
p s for such value of A is summarized in Fig. [5J where 
a cusp-like shape for p s as a function of B emerges: in 
the critical region between 8.51 ± 0.01 < B < 9.25 ± 0.01 
the superfluid phase is present, as testified by the fact 
that here p s is not null. For most values of the mag- 
netic field, modest values of m seem to be sufficient to 
attain good accuracies; close to the border of the critical 
zone, where variations of p s are more sensitive upon in- 
creasing to, higher precision are required though. For all 



the considered values of to, the errors are smaller than 
the symbol size. As an example, for B = 8.6, ranging 
from to = 5 to 40, we obtained values of p s differing 
only by < 5%. By increasing to, indeed we observed 
a vary fast convergence to the asymptotic value of the 
stiffness. This ensured us to obtain reliable results, even 
without pushing further the simulations to larger bond- 
link values. On the other hand, one needs also to increase 
the truncation parameters with to, since too small values 
originate non-monotonic fluctuations in the variational 
energy. — In particular, if an increase of m is not accom- 
panied by a gradual increasing of p and s, the error bar 
in p s increases. 

The scaling behavior of the spin-stiffness is analyzed 
in Fig. [3] for those values of A and B for which there is a 
direct supersolid-to-solid transition. Data are shown for 
the lower critical field at A = 6. Very close to the critical 
field B c the data are described accurately by a power-law 
behavior p s ~ (B — B c )^ a . The value of the exponent 
is very sensitive to the location of the critical point, a 
change in its estimate on the third digit may change the 
value of the fitted exponent up to few percents. By fitting 
all the values up to the vertical bar in Fig. [3] and using 
a value of B c = 8.5052, we get a best fit to the exponent 
of p s — 0.511 which is in very good agreement with the 
theoretical value p s — 0.5 (dashed blue line)^ 

The calculation of the solid order required much larger 
MPS matrix dimensions. However, as already mentioned, 
since for large systems boundary effects are negligible 
when detecting the solid order, we computed S zz (n) by 
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Figure 5: (color online). Spin structure factor S zz (ir) as a 
function of the system size N at A = 5, D = A/2 and for 
different values of the external field. At a value of the field 
B c fH 7.35 ± 0.075 there is an upturn of the curves showing 
that the system becomes solid. 



resorting to a standard variational MPS algorithm with 
OBC, where much larger m values are attainable. To 
guarantee that our data are not qualitatively affected by 
boundary effects, we compared S zz (tt) of Eq. © with 
the one evaluated by summing up only over a fraction of 
the spins corresponding to the central part of the chain 
(say, 1/3 of the total length). The location of the phase 
transition point, where the solid order parameter drops 
from a finite to a vanishing value, do not change, even if 
the value of Osdw inside the solid phase can be different. 

The results for the structure factor are reported in 
Fig. [4] for two emblematic cases. The left panel is ob- 
tained by setting B = 8.8 and A = 6: it corresponds to a 
configuration which is well inside to the cusp of Fig. [5] of 
the supersolid phase. For these values the system should 



hence exhibits a non-null solid order parameter Osdw'- 
this is clearly evident in the left panel of Fig. 21 where 
the value Osdw ~ 4.6 x 10~ 2 is found by extrapolating 
numerical data for N — > oc from the linear behavior in 
N of the quantity S zz (tt). [Notice that the solid order- 
ing can be extracted only for m ~ 100, since at low m 
the data accuracy rapidly deteriorates for larger sizes]. 
On the other hand, the right panel of Fig. 0]is obtained 
for B = 15 and A = 6. It corresponds to a configura- 
tion which is far away from the supersolid region and for 
which the simulations of Ref. 9| predicted that no solid or- 
der should exist (indeed, the system is a supcrfluid there) . 
This is confirmed by our simulations, where we observed 
S zz (n)/N — > in the thermodynamic limit, within nu- 
merical accuracy (Fig. 2J right panel). 

Finally we observe that, for values of the anisotropy 
A < 5.5 in Fig. [TJ there is a direct transition from the 
supersolid to the supcrfluid phase. In this case the tran- 
sition is detected by the vanishing of the solid order pa- 
rameter. In Fig. [S] we show the spin structure factor as 
a function of the system size for different values of the 
external field, fixing A = 5. A scanning of this type for 
different values of A allows to complete the boundaries 
of the supersolid phase. 

In conclusion, we analyzed the supersolid phase in 
a one-dimensional anisotropic spin-1 Heisenberg model 
in a transverse magnetic field, and single-site uniaxial 
anisotropy. By means of an MPS variational calculation 
with PBC, we showed how to determine the spin-stiffness 
and the structure factor, such to locate the supersolid in 
the phase diagram of the system and find the critical 
exponent of the transition to the solid phase. For our 
model of interest, the resulting portion of the phase dia- 
gram containing the supersolid phase is shown in Fig. [T] 
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For large values of A and for D = A/2, the model in 
Eq. |l} at low energies can be mapped onto an effective 
XX spin-1/2 chain in a transverse field, as shown in Ref. 0. 
We analytically extrapolated the critical exponent j3 a for 
such effective model after diagonalizing it in momentum 
space (in presence of a generic twist at the boundary). We 
found a theoretical value $ a = 0.5; this agrees with the 
numerically computed value j3 s , within our accuracy. 



